%Generalized LxF
function u=glxf(u0,dt,maxa,N,al,b,cfl)
u=zeros(1,N);
h=1/N;
la=dt/h;
flux = @(uu) uu.^al/b;
df=@(uu) al*uu.^(al-1)/b;
%Q=0.8
for j=1:N
    if j~=1&&j~=N
        if u0(j+1)~=u0(j)
            Q=cfl^2*(flux(u0(j+1))-flux(u0(j)))^2/(u0(j+1)-u0(j))^2;
        else
            Q=cfl^2*df(u0(j))^2;
        end
        fp=(flux(u0(j+1))+flux(u0(j)))/2 - Q*(u0(j+1)-u0(j))/(2*la);
        
        if u0(j)~=u0(j-1)
            Q=cfl^2*(flux(u0(j))-flux(u0(j-1)))^2/(u0(j)-u0(j-1))^2;
        else
            Q=cfl^2*df(u0(j))^2;
        end
        fn=(flux(u0(j))+flux(u0(j-1)))/2 - Q*(u0(j)-u0(j-1))/(2*la);
    elseif j==1
        if u0(j+1)~=u0(j)
            Q=cfl^2*(flux(u0(j+1))-flux(u0(j)))^2/(u0(j+1)-u0(j))^2;
        else
            Q=cfl^2*df(u0(j))^2;
        end
        fp=(flux(u0(j+1))+flux(u0(j)))/2 - Q*(u0(j+1)-u0(j))/(2*la);
        
        if u0(j)~=u0(N)
            Q=cfl^2*(flux(u0(j))-flux(u0(N)))^2/(u0(j)-u0(N))^2;
        else
            Q=cfl^2*df(u0(j))^2;
        end
        fn=(flux(u0(j))+flux(u0(N)))/2 - Q*(u0(j)-u0(N))/(2*la);
    else
        if u0(1)~=u0(j)
            Q=cfl^2*(flux(u0(1))-flux(u0(j)))^2/(u0(1)-u0(j))^2;
        else
            Q=cfl^2*df(u0(j))^2;
        end
        fp=(flux(u0(1))+flux(u0(j)))/2 - Q*(u0(1)-u0(j))/(2*la);
        
        if u0(j)~=u0(j-1)
            Q=cfl^2*(flux(u0(j))-flux(u0(j-1)))^2/(u0(j)-u0(j-1))^2;
        else
            Q=cfl^2*df(u0(j))^2;
        end
        fn=(flux(u0(j))+flux(u0(j-1)))/2 - Q*(u0(j)-u0(j-1))/(2*la);
    end
    u(j)=u0(j)-la*(fp-fn);
end